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Abstract 

We develop a first-principles electron-transport simulator based on the Lippmann-Schwinger 
(LS) equation within the framework of the real-space finite-difference scheme. In our fully real- 
space based LS (grid LS) method, the ratio expression technique for the scattering wave functions 
and the Green’s function elements of the reference system is employed to avoid numerical collapse. 
Furthermore, we present analytical expressions and/or prominent calculation procedures for the 
retarded Green’s function, which are utilized in the grid LS approach. In order to demonstrate 
the performance of the grid LS method, we simulate the electron-transport properties of the semi¬ 
conductor/oxide interfaces sandwiched between semi-infinite metal electrodes. The results confirm 
that the leakage current through the (001)Si/SiO2 model becomes much larger when the dangling- 
bond (DB) state is induced by a defect in the oxygen layer while that through the (001)Ge/GeO2 
model is insensitive to the DB state. 

PACS numbers: 
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I. INTRODUCTION 


Electron-transport calculations are important tools to investigate and develop materi¬ 
als for new electronic devices. Recently, to obtain more practical knowledge on electron- 
transport properties of nanoscale structures, long-range and large-scale transport simula¬ 
tions have attracted much interest. However, such simulations are very hard task since huge 
computational costs growing with a system size are required. Therefore, it is important to 
develop an efficient electron-transport simulator. 

The Lippmann-Schwinger (LS) equation method proposed by Lang et a/.®^ is one of 
popular methods, which enable us to obtain the scattering wave functions of nanoscale 
structures sandwiched between electrodes by solving the integral equation of the second-kind 
Fredholm equation. When the reference system consists of only bare left and right electrodes 
with the empty transition region, scattering wave functions can be efficiently evaluated for 
a variety of structures of nanoscale junctions set up in the transition region by using the 
same reference Green’s function of the bare electrode system, where the computation of 
the reference Green’s function has only to be performed once. Moreover, for a similar 
reason, the LS equation is utilized in the implementation of self-consistent calculations 
for the convergence of electronic states in inhnitely open systemd®^!. In the conventional 
LS equation method, scattering wave functions are expressed in the Lane representation, 
that is, the LS equation is solved by using a 2-dimensional plane-wave expansion in the 
directions parallel to the electrode surface (lateral directions) and a real-space discretization 
of the coordinate in the direction perpendicular to that (longitudinal direction). In the LS 
equation method, however, one may frequently encounter a numerical difficulty such that a 
part of the Green’s function expressed in a variable-separable form drastically varies due to 
the appearance of evanescent waves exponentially growing and decaying in the longitudinal 
direction. To overcome this issue, in the previous studj^, we proposed the procedure of the 
ratio expression for the Green’s function matrix elements in the Lane representation as a 
remedy for avoiding the numerical collapse. 

So far, we developed the several simulators to elucidate the electronic properties of nanos¬ 
tructures based on the real-space hnite-difference (RSFD) approach!®®, in which the system 
is divided by equally spaced grid points, within the framework of the density functional 
theorjffSI®. For electron-transport simulations, the RSFD method has several advantages 
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compared with the method of the Laue representation from the fundamental and practical 
points of view. Firstly, the hnite differentiation for the kinetic-energy operator is treated on 
the equal footing in all three directions. This avoids numerical errors due to the artihcial 
anisotropy between the lateral and longitudinal directions at any grid spacing. Secondly, 
the computational accuracy can be improved by employing a higher-order hnite-difference 
formula. Thirdly, in the lateral directions, isolated boundary conditions are available as well 
as periodic ones, which enable us to treat electrodes as leads. Furthermore, the algorithm 
of the RSFD method is suitable for massively parallel computin^^. 

In this paper, we present the fully real-space based LS method and the ratio expression 
technique for the Green’s function of the reference system within the approach of the RSFD. 
This method is referred to the grid LS method. To demonstrate the performance of the grid 
LS method, we use it to investigate the electron-transport properties of the (001)Si/SiO2 
and (001)Ge/GeO2 models connected to semi-inhnite electrodes. We also estimate how the 
dangling bond (DB) caused by an oxygen vacancy contributes to leakage currents across 
the interface between the semiconductor and oxide. The results indicate that the leakage 
current attributed to the DB state in the Si/Si02 model is much larger than that in the 
Ge/Ge 02 model. 

In the followings of this paper. Section [IT] gives details of the computational scheme used 
to develop the grid LS method. Section in presents a demonstration of our method, in 
which we use it to examine transport properties of Si/Si 02 and Ge/Ge 02 models and to 
reveal how the leakage current is influenced by the DB state that arises due to an oxygen 


vacancy. Gonclusions are given in Section |IV| and mathematical details are described in 
Appendices A and B. 
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II. COMPUTATIONAL FORMALISM 


We propose an efficient procednre to obtain the solntion of Kohn-Sham eqnation for a 
system where the nanoscale jnnction is sandwiched between semi-inhnite electrodes within 
the framework of the RSFD scheme. The effective potential is close to periodic bnlk poten¬ 
tials as it goes deeply inside the left and right electrodes, so that the whole inhnite system 
can be appropriately divided into three parts: the left electrode, the transition region, and 
the right electrode. The Hamiltonian of the system, H, is dehned by 

if =+ n(r,r'), (1) 

with 


v{r, r') = {vh(r) + v^c(r) + vi(r)}d(r - r') + Vni{r, r'), (2) 

where Vh{r) and v^dr) are the Hartree and exchange-correlation potentials, respectively, and 
vi{r) and Vni{r,r') are local and nonlocal parts of atomic psendopotentials, respectively. 

Assnming that the Hamiltonian in the transition region can be decomposed into an 
nnpertnrbed part and a pertnrbation 6v{r, r') = H — H^, we rewrite the Kohn-Sham 
eqnation as 

/ OO 

dr'6v{r,r')'ijj{r'), (3) 

•OO 

where '0 (t’) is the scattering wave fnnction for an incident wave coming from the left or 
right electrode with the energy E. The snbscript 0 on the variables indicates that they are 
evalnated in the nnpertnrbed reference system. Here, for convenience, we assnmed not 
to contain the nonlocal parts of the psendopotentials. Once the retarded Green’s fnnction 
gd’{r,r'-,E) in the transition region associated with the nnpertnrbed part is known, 
Eq. (|^ is pnt into the LS eqnation in a form of the integral eqnation, i.e., 

^p{r)=^jJ^{r) + jj dr'dr”g^{f’,r'-,E)5v{r\r")ip{'r’") (4) 

with the nnpertnrbed wave fnnction d^{r). Eqnation Q provides a nnihed treatment of the 
Kohn-Sham eqnation and the bonndary condition^^®. In the case where the incident Bloch 
wave 0 *”(r//, 2 :) propagates from deep inside the left electrode, the bonndary condition is 


N 




rj(j)Y^{r//,z) in the left electrode 


N 


( 5 ) 




in the right electrode 


t j=i 
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Here, {r/f, z) is a reflected wave that propagates and decays into the left electrode, 
within the right electrode is a transmitted wave, and rj and tj are nnknown 
reflection and transmission coefficients, respectively. The lateral (x and y) and z directions 
are set to be parallel and perpendicular to the electrode surface, respectively. The system 
is assumed to be periodic in the lateral direction and inhnite in the direction. The case of 
incident electrons coming from the right electrode can be considered in the same manner. 

In this paper, the LS equation is solved within the RSFD schem^^. The RSFD approach 
enables us to treat arbitrary boundary conditions and to calculate the atomic and electronic 
structures to high accuracy. The whole system is composed by the transition region sand¬ 
wiched between semi-inhnite left and right electrodes and is divided by grid points with an 
equal spacing of where and are the length and the number of grid points 

in the y direction (/i = x, y, and z) of the transition region, respectively. Here, we assume a 
2-dimensional periodicity in the lateral directions and employ a generalized z coordinate (k 
instead of Zk, which stands for the group index of z coordinates within the closed interval 
ZkjVf], where TV} is the number of x-y grid planes involved in (k (see Fig. [^; 
A/} corresponds to the order of the hnite-difference approximation for the kinetic energy 
operator in the Kohn-Sham equatiorP^*^ and is chosen so as to include the nonlocal region 
of pseudopotentials to obtain highly accurate results. 

The Kohn-Sham equation is written in a discretized matrix fornP®^ as 

- B{Ck-i; kjjy^iCk-i; kjj) + [E - H{Ck; kJJ)] v[/(a; kjj) - R(a+i; kJJ)^{Ck+^, kjj) = 0, 

(k = -oo,...,-l,0, l,...,cx))(6) 

where H{(k',k^) and B{(k',k^) denoting the A^-dimensional block matrices {N = x 
Ny X TV/) are the diagonal and off-diagonal elements of the Hamiltonian block-tridiagonal 
matrix H{k^), respectively. H{(k',k^) includes the potential on the x-y planes at C = Ck, 
'^{Ck'i is a set of the N values of the wave functions on the x-y planes at C = Cfc, and 
k^ = {k^,ky) is the lateral Bloch wave vector within the hrst Brillouin zone. Hereafter, 
for simplicity, k^ is ignored throughout. We assume that the Hamiltonian in the transition 
region can be decomposed into an unperturbed part H^{Ck) and a perturbation 6V{(k, Ci) as 
well as in the case of the non-discretized treatment mentioned above. When the electrodes 
in the unperturbed reference system are adopted to be exactly those in the perturbed system 
described by H{(k), the perturbation SV{(k,Ci) has nonzero elements only in the transition 
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region (Ci < Cfc(z) < Cm) as 


SViaXi) 


' H{a) - H^ia) 

B{a) - 

Bia.^y - B^a.iY 

0 

V 


{k = i) 

{k = i-i) 
ik = l + l) 
otherwise 


(7) 


Now, by using the discretized retarded Green’s function G^{(k, Ci] E) in the transition region 
associated with the unperturbed part the LS equation is expressed in the discretized 

form as 

m 

= (fc = 0,l,---.m+l). (8) 


which is referred to as the grid LS equation. This discretized form within the framework of 
the RSFD approach unihes Eq. (|^ and the scattering boundary conditions. The boundary 
condition Eq. now reads as 


N 


^(C;t) = 


*^*"(Ca:) + ^G*^j^'^(CA:) in the left electrode {k < 0) 
t=i 


N 


(9) 


t i=i 


in the right electrode (fc > m + 1) 


As the Hamiltonian matrix of the unperturbed reference system, is a block-tridiagonal 
form, the A^-dimensional block matrix E), which is a component of the retarded 

Green’s function matrix G^ = {E — is expressed in terms of the scattering wave 

functions in a variable-separable form as (see Appendix A) 


G^T^CkXuE) 



< D? 


{k < i) 
{k = i) 
{k > i) 


( 10 ) 


Here, {Uy(k)) is the iV-dimensional matrix made of the solutions of the Kohn-Sham 

equation in the case of electrons coming from the right (left) electrode in the reference 
system, that is 


G^Ck) — (^'^%i{.Ck),'^%2iCk), ■ ' ' ) ^?{,Ar(Cfc) j ; 

= (s'IiKO.'I'mKO.--- .S'IkK/.o). 
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( 11 ) 

( 12 ) 








where the iV-dimensional columnar vector (^z,j(Cfc)) denotes the scattering wave 

functions at Cfc for the jth incident wave incoming from deep inside the right 

(left) electrode, where the incident wave is considered to include an evanescent wave as 
well as an ordinary propagating wave; more precisely, taken to be a 

set of the N generalized Bloch states consisting of leftward (rightward) propagating Bloch 
waves and decaying evanescent waves toward the left (right) side, which are the solutions 
of the 2iV-dimensional generalized eigenvalue equatiorP^^I^. The matrix stands for the 
diagonal block-matrix element of the retarded Green’s function matrix, G^{(k,Ck] E), the 
representation of which is derived in Appendix [A| as (see Eq. (AlO)) 


Di = [-B”(a-i)'f/R(a-i)c/H( 4 )-' + ^“(co - B°(a)c'?(a+.)f/ 2 (&)“‘](13) 


with A‘’((Cfc) (—S°((Cfc)) being the diagonal (off-diagonal) block-matrix element of {E — H^). 
Since fo)^(i)(Cfc) includes the exponentially growing or decaying evanescent waves, the 


calculation using Eq. (10) frequently gives rise to the serious numerical error#*. We provide 


a remedy for this problem as follows. 

Introducing the ratio matrices and Y)? at two successive (k points, which are defined 


as 


W° = C/h(4-i)(C/r(C0)’‘ (14) 

n" = i7»(&+.)(c/j(a))"‘. (15) 


respectively, we obtain the following (m -|- 2)-dimensional block-matrix expression for the 


retarded Green’s function Eq. (10): 


G "' rO 
T 


m+1 


no 

-^0 




0 

1 

IlT-DS 





m+1 


i=i 

m +1 


X»£>» 


n-T-'Ds.+i 

J=2 


no 

^m+l 


]=m 


j=m 


j=m 


( 16 ) 









that is, we rewrite the block-matrix element of in Eq. (10) as 


( I 

(*: < i) 

j=k+l 


GJr^{CkXi;E) = { Df 


{k = l) , 


{k>l) 

K, j=k-l 


and from Eqs. (13)~(15) the diagonal block-matrix element reads as 


Dl = + A«{Ck) - B“(4)n°] 


01 -1 


( 17 ) 


( 18 ) 


In the following subsections |II A| and |IIB[ we will give efficient numerical calculation tech¬ 
niques for the ratio matrices {X^} and {Y^} without employing the matrices {h^/i;(Cfc)} 
and {Ul{(k)} which include evanescent waves explicitly. Our previous studj^ verihed that 
the introduction of the ratio expression such as Eqs. Q-Q into the retarded Green’s 
function enables us to avoid the numerical collapse originated from the appearance of the 
rapidly growing and decaying evanescent waves. By contrast, in LS simulations of elec¬ 
tron transport through long conductor systems using the conventional Green’s function in a 
variable-separable form, the numerical collapse is inevitable. 

In the solving of Eq. ([^ by using the iterative method such as the conjugate gradient 
method, the operation of X (Cfc) 0')d'(C/0 in Eq. (|^ is carried out as follows: 

m 

Y. G?(4. Q)SV{(,t, CO'S (Cl') = 'S'(Ck) + PdCt) + Pr{C0 

l,V=l 

{k = ,m + l), (19) 


where 

m 

( 20 ) 

1=1 

n(co= ny"'>i''(Co)+ nK/.>i<'(Ci) + --- + yyi>i''(a-i). (21) 


pock) =w\,>i>'(4+i)+w^,A^,<i<'(a+2)+■ ■ ■ + n y“' ^'(Cm+i). ( 22 ) 

j=k+l 

It is easily shown that the sequences {PLiCk)} and {PniCk)} satisfy the following recursive 
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relations: 


PdCk) 

PniCk) 


0 

(A; = 0,1) 

(23) 

i)?_.|Pi(a-i) + >t'(4-i)i 

(A: = 2,3, • • • ,m -t-1) 

0 

.’f^.|Pfl(&+i) + 4''(a+.)] 

[k = m + 1, m) 

(A: = m — 1, • • • ,1,0) 

(24) 


Here, we used the fact that SV{Ck,Ci) = 0 outside the transition region of < Ck{i) < Cm- 
It should be emphasized that since the elements of in Eq. (17) are no longer in a 
variable-separable form, the amount of [{m + 2)N]‘^ for each multiplication is expected to 
be required; nevertheless, it is reduced to the order of (m -|- 2)iV^ by virtue of Eqs. (23) 


and (24), which means that the present method does not suffer from the numerical collapse 
without increasing the computational cost. 


A. Jellium Electrodes 


The case in which electrodes are approximated by structureless jellium models is treated. 
The jellium electrode approximation has been successfully applied to the interpretation of 
electron-transport properties with less computational loacP®®. A free electron system is 
chosen as the unperturbed one with the Hamiltonian where a completely flat potential 
is assumed, for simplicity. The Green’s function in the free-electron system is more con¬ 
veniently described by using Zk instead of (k- In Appendix we discuss the analytical 
expression of the Green’s function in terms of Zk in a general A// case. 

We here give details on the implementation of the analytically expressed retarded Green’s 
function in the 3-dimensional central hnite-difference (7V/=1) case, for example, which is 
written by 


G't (G/j’ r//j>, zi; E) 


— ^ ^ / R \ 

) " ^Ili3 '^11 ■>3^') ^l\) 

iN ^ ^ sin/Ci/i^ 

(25) 


Here, 


G//,n — {Gn^,Gny) — 


27r 27r 

hyNy 


n, 


(26) 
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^//j ~ i^jx^Vjy) the lateral coordinates with jx{y) = 1,2,-•• ,N^{y) [Nx{y) is chosen an 
odd integer for convenience], and 
1 


cos 


-1 


/Cl = <j cosh-^ (l -hl{E- 

TT + i cosh"^ T-l + hl[E - 




■■■ ^ 


E <E, 


(A/'/ = l) 


(27) 


V 


1) ^ 

' ' ' "r ^2 — 


with 


^ {1 - cos + /cf) /l^} + {1 - cos {Gn^ + /cj) /zj , 




(28) 


In the derivation of Eqs. (25)-(28), we used Eqs. (B22) and (B23) and the extension of 


Eq. (B27) to the case of the 3-diniensional space. 

Since dehned by Eq. (17) is the diagonal block-matrix element of the retarded Green’s 
function G^{zk, zi; E), the jth row and j'th column element {Dl)jj/ is expressed as 

= G^rir/fj, Zk, Zk] E) 

1 


K 

iN 


Nx-1 

2 2 2 


exp(z(G//,„ + kjj) ■ {r//j - r//j,)) x ^ 


JV!L=1„ _ Ny-1 
2 "^y~ -2— 


sin/Cl 


(29) 


One can see from Eq. (29) that {Dl)jj/, and thus D^, is /c-independent owing to the trans¬ 


lation invariance in the z direction. On the other hand, by Eq. (17), and are given 
by 




(30) 

(31) 


and from Eq. (25), the jth row and j'th column matrix element of G‘!^{zk±i, Zk] E) are 
described as 


(G5?(zfc±i, Zk, E)). ., = (r//j, Zk±i,r//j,, Zk, E) 

Nz,-1 ^H-1 

= ^ i i exp(z(G//,„ + fc/f) 

Nx-1 Ny-1 


m/.i - '•he 


)) 


.V, _ JVa; —1 iV^- 

rix- -^ ny = - \ 


sin JCih. 


exp(z/Ci/z^). 


( 32 ) 


This implies that X)! and are also /c-independent and 




( 33 ) 
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After some calculations, 




Ny-1 


N 


Y1 + ^//) ■ “ ^//J')) X exp{ilCih,) (34) 


Nx-1 Ny-1 

ny= -L 


is obtainecP^. Hereafter, and are denoted by and X°, respectively, since they are 
A;-independent. 

The products of the matrix (-A°) and vectors {f{r//j, Zk)\ j = 1,2, ■ ■ ■ , as required 


in the computations of Eqs. (20), (^23^ and (24) can be easily carried out in the momentum 


space, since they are written in the convolution form of the 2-dimensional discrete Fourier 
transform. Owing to the orthogonality of the plane waves, the Fourier transformed and 
are represented as the diagonalized matrices, i.e.. 


^ = exp(-z(G//,„ + kjj) ■ r//,, + + kfi) ■ r//,y)(D°)y,- 

1 


^//,i ^//,y 


— ih^Snn' ■ . / T \ ) 
sm(/Cih^) 


^ n.n' + fe//) ■ r,„ + + kf) . 


A/j A/j' 


= inn' ewiXih^), 


(36) 


(36) 


respectively. Finally, one can obtain the matrix elements of the Fourier transform of the 


terms shown in Fqs. (20), (23) and (24) as 


/ 

0 


y^6V{zk,zi)^{zi) 






(37) 

(A) = 0,1) 

{k = 2,3,--- ,m+l) ’ 

(38) 

{k = m+l, m) 

{k = m — l, ••■,1,0) ’ 

(39) 


respectively. 
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For calculating the product in Eq. ([^, the computational cost of 0{Nin x x iV^) is 
required. However, by introducing the 2-dimensional discrete fast Fourier transform (FFT) 


algorithm, the cost of the product shown in Eqs. (37)~(39) decreases to 0(iVj„ x x 
N\ogN), since the off-diagonal elements of the Fourier transformed matrices J^[D^] and 


[X°] are zero as seen in Eqs. (35) and (36). The Fourier transform of a columnar vector 
and the inverse Fourier transform of [D^] x [fk] and [X°] x X [fk] are carried out at 
each Zk point using FFT algorithm. Here, X [fk] represents the Fourier transformed vector 
of [fir//.. j) ^k)\ j = 1, 2, • • • , X}. Thus, the maximum order of the calculations is improved 
from 0{Nin X X| x X^) to 0{Nin x x XlogX). The above mentioned discussion on the 
central hnite-difference approximation can be straightforwardly extended to the cases of the 
higher-order hnite-difference approach. 


B. Crystalline Electrodes 


A general case is discussed where a system with atomistic crystalline electrodes is chosen 
as the unperturbed reference system; one electrode is confronted with the other across the 
empty transition region. We present efficient procedures for calculating the ratio matrices 
X° and in this case. 


The matrices Xq and dehned by Eqs. (14) and (15) are described as 



= (^"(C-OT'E'lKo). 

(40) 



(41) 


where X]l°(Co) (X]fl(Cm+i)) is the self-energy term dehned on the left- (right-)electrode 
surface and can be calculated by using the continued-fraction equation; for the details of the 


derivation of Eqs. (40) and (41) and the computation of the self-energy terms, see Refs. 10 
and 18. For the sake of comparison, we note that UR(L){Ck) dehned by Eqs. 0 and ( p!^ 
is identical to Q'^^^\Ck) of Eq. (15) in Ref. 18. We also emphasize that the accuracy of 
'^R(L){Ck) is enhanced by making use of the continued-fraction equation in a self-consistent 
manner, as shown in Eqs. (16)-(18) in Ref. 18. 
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It should be noticed that the terms {X°} can be sequentially computed as 


A'? = (A“(Co) - -S“(C-i)'a;) ■'«“((()) 
X'i = (/l»(Ci) - B“(Co)'A'!’)-‘b»(Ci) 
As" = (yl"(C2) - B"(Ci)'A?)"‘b"(C2) 


A2.+, = (2l»(C„) - B»(C„-,)tA»)-'B»(C„), (42) 


which are easily derived from Eq. (A3), and similarly, the iterative series of {1^°} are ob¬ 


tainable from Eq. (A4) as 


r^l = (A"(C™) - 

r™-2 = (A"(C™-.) - B"(C™-.)n",_,)“'B"(c„-2)' 


= (A"(Ci) - B“(Ci)F»)-‘B»(Co)t. 


(43) 


The recursive relations Eqs. (42) and (43) allow us to calculate all the matrix elements 


by a linear scaling operation (order-A^ calculation procedure) at a limited computational 


cost. It is also noted that using Eqs. (42) and (43), and are stably computed without 


involving error accumulation since the errors due to the appearance of evanescent waves are 
eliminated by introducing the ratios of these waves at two successive grid points. Finally, 


the diagonal block-matrix element is given by Eq. (18). Once D^, and Y)P for any k 


{0 < k < m + 1) are known, all of the matrix elements of in Eq. (17) are determined. 


and the algorithm of Eqs. (23) and (24) can be utilized. 
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III. APPLICATIONS 


To demonstrate the performance of the grid LS method, we examine the electron- 
transport properties of models of semiconductor/insulator interfaces sandwiched between 
semi-inhnite electrodes. Recently, the germanium-based metal-oxide-semiconductor field- 
effect transistor has attracted signihcant attention because the electronic band gap of 
germanium (~ 0.66 eV) is lower than that of silicon (~ 1.12 eV), which allows for reduced 
operating voltages. In a highly integrated circuit, it is known that a large leakage current 
is induced by defects such as impurities and oxygen vacancies in the thin gate oxide layer. 
So, the relationship between the DB introduced by defects and leakage current in Si/Si02 
interfaces has been extensively investigatecP^®, while the role of the DB state in Ge/Ge02 
interfaces is controversial. One of the present authors (T. O.) has performed several in¬ 
vestigations on Ge/Ge02 interfaceP®^. In recent work, the relationship between atomic 
conhgurations and electronic structures of (001)Si/SiO2 and (001)Ge/GeO2 models with 
DBs was explored using first-principles simulations within the framework of the local den¬ 
sity approximation (LDA)3S1. It was found that the Si-DB state is located near the midgap 
of the Si substrate corresponding to the Fermi level, while the Ge-DB state lies near the top 
of the valence band which is 0.3 eV below the Fermi leveP^. 

To examine how DB states with different characteristics affect leakage currents, we per¬ 
formed transport simulations of electrons flowing across the (001)Si/SiO2 and (001)Ge/GeO2 
models. The magnitude of the leakage current flowing through insulators is so small that 
it can be easily affected by interactions between electrodes and interface models and by the 
value of the energy band gap, which is underestimated by the LDA calculation. Therefore, 
in this paper, we discuss qualitatively the ratio of the leakage current between models with 
and without a defect. 

Figure illustrates a unit cell of each interface model. In these models, the side lengths 
of the cell in the lateral direction parallel to the interface for the Si/Si 02 (Ge/Ge 02 ) model 
were taken to be the experimental lattice constant of bulk Si (Ge), Oq = 5.43 (5.65) A. The 
thicknesses of the Si02 (Ge02) layer and the Si (Ge) substrate were 7.34 (7.25) and 7.18 
(7.37) A, respectively. In calculations for the models with an oxygen vacancy, we introduced 
the defect into a supercell comprising 4x4 unit cells in the lateral direction (Fig. [^; this 
is large enough to avoid interactions between defects in neighboring cells. Two Si (Ge) 
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DBs were generated near the interface between the Si (Ge) substrate and the oxide layer in 
one unit cell by removing a bridging oxygen atom in a manner similar to that used in the 
previous studj^. One of the DBs is passivated by a hydrogen atom, while the other remains 
with the Si (Ge) atom of the center back-bonded to two neighboring Si (Ge) atoms and 
an oxygen atom (•Si=Si20 (•Ge=Ge20)). For the no-defect models, the unit cell of each 
model, depicted in Fig. was employed with 4x4 sampling fc-points in the 2-dimensional 
Brillouin zone for comparison with the models having defects. 

We hrst optimized the atomic and electronic structures of the models. First-principles cal¬ 
culations based on the RSFD approach were performed in the manner described in Ref. [42] 
with a grid spacing of 0.15 A. The size of the supercell in the [001] {z) direction was taken 
to be 5ao, including a large enough vacuum region, and the top and bottom layers of the 
models were terminated by hydrogen atoms. As shown in Fig. which illustrates the re¬ 
laxed conhgurations, the Si atom with the DB is pulled down to the Si substrate whereas 
the Ge atom with the DB is slightly raised toward the oxide layer. 

Next, we examined the leakage currents caused by introducing the oxygen vacancies into 
the models. Employing the optimized effective Kohn-Sham potential, we used the grid LS 
method to evaluate the scattering wave functions for electrons incident from the bottom-side 
electrode. The conductance at the limits of zero temperature and bias are described by the 
Landauer-Biittiker formulaP^. In the transport calculation, the top and bottom sides of each 
model were connected to aluminum jellium electrodes without terminating hydrogen atoms. 
The Wigner-Seitz radius was 2.07, which corresponds to the valence electron density of 
bulk aluminum. 

Figurej^shows the computed conductance spectra for the no-defect Si/Si02 and Ge/Ge02 
models as functions of incident electron energy measured from the valence band maximum 
(VBM) of the substrate. Although some small peaks derived from the bulk states in the 
valence band of the Si and Ge substrates appear in Fig. both models exhibit highly 
suppressed conductivities in the band-gap region between the VBM and conduction band 
minimum (GBM) of the substrates. Figure represents the conductance spectra for the 
Si/Si02 and Ge/Ge02 models with the oxygen vacancy. No remarkable peaks appear in 
the spectrum of the Ge/Ge 02 model; however, for the Si/Si 02 model, a peak with high 
transmission occurs around VBM -|- 0.41 eV, where electrons flow through the oxide layer 
via the Si-DB state as shown in the charge density distribution of the scattering electron 
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(Fig. [^. In addition, Fig. [^exhibits contour plots of the density of states (DOS) integrated 
in the lateral directions (left panels) and the charge density distributions of the DB states 
(right panels) for the defect-introduced models without electrodes. In Fig. [^a), the white 
arrow identihes the peak in the DOS derived from the DB state between the VBM and CBM 
of the Si substrate. In Fig.j^b), the Ge-DB state is coupled to the states in the valence band 
of the Ge substrate. The relative positions of the VBM, CBM, and DB state are modulated 
when each model is connected to electrodes. When a DB appears in diamond-structured 
semiconductors, there are two possibilities: the DB state tends to become either more s 
type or more p typ^^. In the p type DB state, the three remaining bonds tend to become 
sp‘^ hybridized and, to reduce the strain, prefer to be in a plane. This occurs in Fig. [^a) 
wherein the Si atom with the DB is pulled down and the DB state spatially extends in the 
[001] (z) direction. This behavior degrades the insulating properties of the Si/Si02 model. 
In contrast, when the DB state is inclined to be an s type, the three remaining bonds tend 
to become p types. In this case, the angular separation of these bonds is reduced from that 
in the tetrahedral structure where the separation angle is 109.5°. As a result, the atom with 
the DB moves away from the three bonded atoms. Therefore, for the DB state of the raised 
Ge atom, the charge density of the state is distributed in the lateral directions compared 
with that of the Si-DB state, and the Ge-DB state is coupled with interface states of the 
Ge substrate (Fig. [^b)). This behavior barely contributes to electron transport across the 
model. Gonsequently, by introducing the oxygen vacancy, the leakage current in the Si/Si02 
model increases by a factor of 162.9, while that in the Ge/Ge02 model increases by a factor 
of 11.^. 


IV. CONCLUSION 

We have presented the grid LS equation method based on the fully real-space algorithms 
to elucidate the scattering wave functions in nanoscale structures sandwiched between semi- 
inhnite electrodes. It is shown that the numerical collapse due to the exponentially growing 
and decaying evanescent waves and the computational costs can be restrained by using the 
ratio expression of the retarded Green’s function obtained analytically (jellium electrode 
case) and by incorporating the self-energy matrices and applying the recursive formulas 
to the ratio matrices (crystalline electrode case). To demonstrate the performance of our 
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method, we used it to calculate the transport properties of (001)Si/Si02 and (001)Ge/Ge02 
models attached to semi-inhnite electrodes. The results show that the DB state in the 
Ge/Ge 02 model gives a much smaller contribution to leakage current than that in the 
Si/Si02 model. Our procedure can precisely and efficiently extend knowledge of the physics 
underlying the transport of electrons through nanoscale structures. 
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Appendix A: Variable-separable-formed retarded Green’s function in the RSFD 


approach 


In this appendix, the subscript 0 denoting the reference system is omitted, for simplicity. 
Let us consider the product of the matrix {E — H) and the /th columnar vector of G^{E), 
Cd -S)} {k = ■■■ ,l — 1,1,1 + 1, ■■ ■). The retarded Green’s function is constructed 
from outwardly propagating and decreasing waves, and then, taking it into account, we 
assume this columnar vector to be represented by 

• ■ , Un{Ci-2)Gn{Ci), Ur{Ci.^)Gr{Ci), D{Ci), Ul{Ci+i)Gl{Ci), t^L(0+2)G^(0), •••]*, (Al) 

where Gr(l){Ci) and D((i) (= G^{(i,(i; E)) are unknown block matrices, and Ur(l){(i) is 


defined by Eqs. (11) and (12). Hereafter, Ur(^l){Ci)i Cr(l){Ci) and D{(i) are abbreviated to 


^ ^R(L) respectively. By definition, the abovementioned product satisfies 
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— Bi-2 
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Ul^-2Cl^ 
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-bU 

Ai-i 

-Bi., 







0 



-bL 

Al 

-Bi 




Di 
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I 




-B\ 

Ai+i 

—Bi+i 



Ut+iCt 
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0 



-d+i 

Ai+2 



Ut,2Gt 


0 

V 








\ 



V/ 


the /th, (A2) 


where Ai = E — H{(i) and Since is a set of the solutions of the 

Kohn-Sham equation, the following equations hold: 


^kUk+l — 0, 
BkUk+i = 0, 


bUuI, + 

-d-it'hi + Awi 

(k = ■ ■ ■ ,l — 1, /, / + 1, • ■ ■) 


(A3) 

(A4) 


From Eqs. (A2)-(A4), one sees that the unknown matrices and Di are required to 

satisfy the equations 


- Bi.^Di = 0 , 

+ A,D, - = I, 

-BjD, + BjUl-Cl- = 0, 
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(AS) 

(A6) 

(A7) 














R( T') 

and thus, Eqs. (A5) and (A7) lead to the relationships between C, ^ and Di as 


KV'di, 

(AS) 


(A9) 

+ A, - BiUt+Myk■ 

(AlO) 


and Eq. (A6) decides Di to be 
Di = 

In consequence, the retarded Green’s function Cil E) can be described in the following 

separable form: 

' U^{U^)-^Di {k<l) 

Di 


G^CkXi;E)={ 


{k = i) 


(All) 


utiuty^Di {k>i) 


Appendix B: Analytical expression of Green’s function for free electron system in 
the RSFD approach 

A 1-dimensional system is firstly considered for simplicity. In the RSFD approach, the 
kinetic-energy operator K = — is represented by the matrix and the kinetic- 

energy term in the Kohn-Sham equation is written as 

+ G_j^j+i'iIj{zi_m^+i) + ■ ■ ■ 

■ ■ ■ + + GqiIj^zi) + Gi'0(^£+i) -!-■■■ 

■ ■ ■ + CV/-it/’(2;f+A/'/-i) + CV/t/’(2;r+A/>) ) (Bl) 

where A/} is the order of the finite-difference approximation, h is a grid spacing and the 
weight coefficients Gi {i = i — J\ff, i — Mf -|- 1, • • • ,i + Mf) are determined using the Taylor 
expansiorP^. 

In the A//th order finite-difference approximation, the Green’s function matrix is 

determined as satisfying 

(^Z - G^^f\Z) = J, (B2) 

where Z is a complex number and I is the unit matrix. The £ th row-£' th column element 
of the Green’s function matrix G^-^^^Zk, zi] Z) is described in a spectral representation as 

(B3) 


'7^ 4>p{^k)4>p{zi) 

G'-^ f>{zk,zi]Z) = I — — zzjTjydp, 


' — n/h Z — E- 
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where and 4>p{zi) are the eigenvalue and eigenvector of \ respectively, obtained 

by solving the eigenvalue equation (ppi^zi)^ and are given by 


E: 


{Mf) _ 


Nf 


2h2 


Co + 2 Cm COS mph 


m=l 


(B4) 




M~-l) = 7^ 



^ipzi 


27r 


where — | < p < and (j)p{ze) is normalized, i.e., 

OO 

E, 4‘‘r(Z]4>r'(^t) =S{p- p')- 


(B5) 


•^=—00 


Substituting Eq. (B4) into Eq. (B3) and changing the integration variable from p to 9 = ph 


and subsequently from 6^ to a; = e®^, we obtain 




A/> 


-d9 


K^Z + -Co + Cm COS mO 


2Tii 


m=l 




-du. 


f 


(B6) 


h^Z + C„ + t ^ C„(l7"' + w-’”) 


m=l 


The integration can be carried out along the unit circle in the complex u plane based on the 
residue theorem. In the following, we introduce a sensible manner of picking up the poles 
inside the unit circle that contribute to the integration. These poles ca’s are the A// solutions 
of the equation 


1 1 

h^Z + Co+ -Y, = 0- 


m=l 


We now define a new variable s as 


and rewrite Eq. (B7) as 


s =-{u + u ^) = cos6, 


h'^Z + -Co + Cm cos mO = 0, 

m=l 


(B7) 


(B8) 


(B9) 


which is the A//th order algebraic equation with respect to s and its solutions are denoted 
by Sn (n = 1, 2, • • ■ , A/}). Resultantly, for each s„, the poles cn^’s are given by the solutions 
of the quadratic equation — 2snUJn + 1 = 0 as 

= Sn ± (BIO) 
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When the imaginary part of Z is non-zero, either of or cai ^ is inside the nnit circle, 
while the other is ontside it. Hereafter, the inside pole is defined as a;„, and the other is 


obtained by Using the residne theorem, the integration of Eq. (B6) is carried ont to 
yield 


2/,2 




Mj 


{u^n (^n ) 


m = l 
m^n 




K 


f 




OJn 


^ at, 

(^n ^m) 

m=l 

mj^n 


(Bll) 


Finally, introducing /C„ by s„ = cos/C„h, hence, Un = e*^"^ (Im{/C„} > 0), we obtain 


G^^f\zk,zi-Z) 




Nf 

E 


I I 


f n=l 


Nf 


sin JCnh (cos ICnh — cos /Cm/l) 


m=l 

m^n 


(B12) 


The retarded Green’s function is given by 


zi-E) = hm G^^i\zu, zu E + ie). (B13) 

£->0 + 

It is noted that /C„ is a multivalued function since it is defined a.s ICn = \ cos“^ s„. The 
branch of cos“^ should be chosen to satisfy the requirement of Im{/C„} > 0, which guarantees 
that Un exists inside the unit circle in the complex u plane. Thus, the following relationship 
is established: 


f If Im{s„} < 0, then 0 < Re{/C„} < — 

< ^ . (B14) 

If Im{s„} > 0, then < Re{/C„} < 0 

Proof: Consider complex numbers s = ^ + ir] and K, = k + in with the relationship of 
s = cosJCh in the interval of \k\ < It is readily shown that when rj < 0 and k > 0, there 
exists the one-to-one correspondence between {s} and {JC}. In this case, /C is uniquely so 
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defined as to satisfy Im{/C} = k > 0 and the relationship 

1 


cos ^ 7 


e>o 


k = 


h 

^ ^ cos“^ (-7~) ■ ■ ■ ^ < 0 , 


K, = — cosh ^ 7 ’*' 
h 


where 7 ^ is dehned as 


j+ 7^ + 1 ± \/(^^ + yf - 1)^ + 4^2 


(B15) 


(B16) 


and cos ^ (cosh in Eq. (B15) is the principal value of the inverse trigonometric (hyperbolic) 


cosine function. Thus, k varies as 


/' , 71 

0 <k< — ■■■^>0 
- 2h ^ - 




(B17) 


In the derivation of Eqs. (B15)-( B17| ), we used well-known formulas, 

cos {k -f iK)h = cos fchcosh nh — i sin fchsinh nh 
sin^ kh + cos^ kh = 1 ■ (B18) 

— sinh^ nh + cosh^ nh = 1 

On the other hand, in the case of 7 > 0 and k > 0, then k and k are determined in the 
same manner as 

/ 1 . 

... e>o 


1 -1 - 

— — cos 7 
h 


k = 


-^cos ^ (-7 ) ... ^ < 0 , 


(B19) 


K = — cosh ^ 7 ’'“ 
h 


respectively, and k varies as 


f tt 

-- < fc < 0 

2h - 


71 , 71 

-- < k < -- 

h 2h 


e>o 

e<o 


(B20) 
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(Q.E.D) 

It is straightforward to extend the above argument to the 3-dimensional case. We deal 
with a free electron system in which the discretized space is infinite in the z direction and 
periodic in the x and y directions. The Green’s function in this case is described in a spectral 
representation by 

Nx — 1 ^ ^ 

4’nx,ny,p{'f'//,j: Zk) 0 n 2 :,n„,p(^//j'5 


ir//J: Zk, Vn^y, zi;Z)= 


^ 2 ^y— 2 


Z -E.. 


{Ns) 

rix.ny.p 


-dp, 


(B21) 


where Eli^fny,p and (j)n^^ny,p{'f’//,j, Zi) are the eigenvalue and eigenvector of the 3-dimensional 
kinetic-energy matrix, respectively, i.e.. 


pW/) _ pW/) pW/) 

^Tlx,ny,p — ^nx,ny “T 


E, 


{Nf) 


Nf 


2hl 


2hl 


E, 


i^f) ^_L 

2hl 

<Pnx,riy,pit'll,j, Zi) = 


Co + 2'^ Cm COS mCnx hx 

m=l 

Nf 

([ 7 q “h 2 ^ ^ Cyii cos ttiG yi^hy 

m=l 
Nf 

C{) + 2 Cm cos mphz 

m=l 

exp (iG//^n ■ r//j + ipz^ 


(B22) 


2-KNy,Ny 


Here, = ixj^,yjy) are the lateral coordinates with jx(^y) = 1, 2, • • • , W(y) [for convenience. 


Nx[y) is chosen an odd integer], and the dehnition of is same as shown in Eq. (26). 
Now, the Green’s functio] 

G^^^\'f'll,j^Zk,riij,,zi-,Z) 


Now, the Green’s function represented by Eq. (B12) reads as 

hi 


i2^f-^N^NyCMs 

iVx-l 
2 


X 


i ■ irii,j - r//,y) + I Zfc - Zz I) 


ri=l 

"■X 2 “H— 2 


Nf 


sin ICnhz (cos ICnh^ — cos ICmhz 


m = l 
m^n 


(B23) 


where /C„ = ^ cos ^ under the requirement of Im{/C„} > 0, and is the solution of the 
Affth order algebraic equation with respect to s (= cos 6 ') 

1 Nf 


hl(z - EiCfl) + N + Y1 Cmcosrne = 0 . 


(B24) 


771=1 
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In the following, we present the analytic representation of the retarded Green’s functions 
in Aff - I ~ 4 cases in the 1-dimensional free electron system; there exists no analytic one 


in a case of A// > 5, since the algebraic equation (B9) with A/} > 5 can not be solvable 


analytically according to Galois theory. Given the solutions of Eq. (B9), irjn {n = 

1 , 2, ■ ■ ■ , A//), Kn = kn+iHn are determined from Eqs. (B14)-(B20), and hnally we obtain the 


analytic representation of the Green’s function Eq. (B12). Hereafter, we choose Z = E+ie {e: 


an inhnitesimal positive number) in Eq. (B12) so as to deal with the retarded Green’s 
function. 


1. case of central finite difference (A// = 1) 


Substituting Cq = —2 and Gi = 1 into Eq. ( |B9[ ), we have the equation 

s — 1 -I- h‘^{E + is) = 0, 

and its solution 

Si = = 1 — K^E — ih'^e. 


(B25) 


(B26) 


Since Eq. (B26) indicates that r/i —)■ 0 (an inhnitesimal negative number) in the limit of 


e —)■ O'*", ICi = ki + in is determined from Eqs. (B15) and (B16) such that 


fci = — cos ^ (1 — h‘^E ) , Ki = 0 


ki = 0 
u ^ 


, Ki = — cosh ^(l — h^E) 

, Ki = Y cosh“^(—1 -1- h‘^E) 
h ^ ' 


^<E<- 
E < 0 


(B27) 


2. case of 5-point finite difference (A// = 2) 


Substituting Co = —5/2, Ci = 4/3 and C 2 = —1/12 into Eq. (B9) yields the quadratic 
equation with respect to s, 

s^ — 8s -I- 7 — 6h'^{E -I- ie) = 0. 

This equation has two solutions si and S 2 given by 
Si = + ir]i = A - y/9 + 6h‘^{E + ie), 

^1 = 4 - V 9 + Qh^E , 771 ^ 0 - 

.^1 = 4 , Pi = — -y]9”+l3h2^ < Q ... ^ < 


(B28) 


(B29) 


where 


2h2 


< E 
3 ’ 


2h2 
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and 


where 


S 2 = 6 + *h 2 = 4 + ^/9 + 6h'^{E + ie), 
6 = 4 + \/9 + Qh'^E , 7^2 0 + 


(B30) 


2/i2 


< E 


6 = 4 


, ri2 = \/\9 + 6h'^E\ >0 • ■ ■ E < 


2h2 


Subsequently, from Eqs. (B15), (B16) and (B19), /C„ = kn + ii^n can be described in an 
analytic form. 


3. case of 7-point finite difference (A// = 3) 


The substitution of Co = —49/18, Ci = 3/2, C2 = —3/20 and C3 = 1/90 into Eq. (B9) 
leads to the cubic equation with respect to s. 


o 27 2 109 45,2,^ . , ^ 

s -—s + 33s--—I— —h [E + ie) — 0, 

4 4 2 


whose solutions are determined according to Cardano’s formula as 


(B31) 


Si = 6 + where — + t and r/i —)■ 0 , 


+ T 


T ' — T . , . 

S2 = 6 + where 6 =-^- and r]2 = - - - > 0, 

T+ — 7-“ 7-+ -|- 7-“ 

■S3 = 6 + where 6 =-- and 773 =---< 0. 


(B32) 

(B33) 

(B34) 


Here, 


4 


V(144h2E + 155)2 + 5 X 193 ± (144h2E + 155)|, (B35) 


and JCn = kn + iHn can be analytically represented using Eqs. (B15), (B16) and (B19). 


4. case of 9-point finite difference (A// = 4) 


Now, Substituting Cq = —205/72, Ci = 8/5, C2 = —1/5, C3 = 8/315 and C4 = —1/560 


into Eq. (B9), we obtain the quartic equation with respect to s. 


_ ^^3 ^ 27s2 -^s+ — - 70h\E + te) = 0. 
9 3 9 ^^ 


(B36) 
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Ferrari’s solutions to Eq. (B36) are utilized. After tedious but straightforward calculations, 
we have the following representations: 


. 16 ^ (3 

6 = y + cr - Y a ’ 

. 16 


, r;i--Wc7^ + a- l<0 


a 


^ 16 / 9 n + 

6 = vr + cr + \/ -cr - a + - , -t 0+ 
9 V (T 

. 16 

6= g-+a 


, j/2 — \/ + a - ^ > 0 


(B37) 
' ■ Eq < E 

■ ■ E <Eo 

(B38) 

Eq < E 
E <Eq 

(B39) 


16 / 9 n 

^3 = — - a , r]Q = + a + - 

9 V cr 

. 16 r. ^ 

^3 = y- a , r]Q = \^a^ + a + ->Q 


Here, 


. 16 / , /3 

U = — - a , 774 = \/o-^ + a + ->0 

9 V cr 

. 16 rr ^ 

44 = — - cr , ?74 = -Wcr^ + cr + - < 0 
9 V cr 


7-31 , 23-7-181 

cr = ::—::7r , p = 


E' <E 


E<E' 


E' <E 


E<E' 


(B40) 


2 ■ 33 


36 


cr = < 


cr 


cr 


-^ + \JP+ \/i^ + - \J-p + \/ic 3 + 


+ \/p + + p^ + Y^p - \/z^^ + 

o 


0 < z/ 

z/ < 0 


(B41) 


z/ = ao + aih‘^E , p = -z/ + 02 , 


En = 


Oq — 


Z/q — Uq 
aiK^ 

32 ■ 5■19 


, E'q = 


^0 - Qq 
aiK^ 


, Oi — 


2 ■ 33 ■ 5 


5 CX 2 — 


32 ■ 5 ■ 3623 


22.312 ’ 7.312 ’ 2 - 7-313 ’ 

and vq is the solution of = 0, which is evaluated as z/q = — 0.3563, and Vq = 

— ^(03 + 1/3)3 — 202^/303 = — 0.4132 with 03 being l/2(P/2a. The usage of Eqs. (B12)- 
( B20[ ) leads to the analytic representations of /C^ = + iKn and the retarded Green’s 

function. 
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The above treatment for analytically describing the retarded Green’s function is readily 


extendable to the 3-dimensional case using Eqs. (B22)-(B24). 
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In each model, the leakage current was calculated by integrating the tunneling current over the 
energy ranges between VBM and CBM of the substrate. 
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FIG. 1. Sketch of the relationship between z and ( in the computational model in the case of JVf = 
3. The system consists of the transition region sandwiched between the left and right semi-infinite 
crystalline electrodes. In the left electrode, the incident wave and the reflected waves consisting of 
the propagating and evanescent ones are illustrated by <!>*"■ (Cfc) and respectively, while in 

the right electrode, the transmitted waves composed of the propagating and decaying evanescent 
ones toward the right side are denoted by Here, x(y) and z are coordinates perpendicular 

and parallel to the nanoscale junction, respectively. 



FIG. 2. (Golor online) Schematic views of unit cells of (a) Si/Si02 and (b) Ge/Ge02 models. 
Dashed lines indicate the boundaries of the cell. White, blue, green, and red spheres represent H, 
Si, Ge, and O atoms, respectively. 
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FIG. 3. (Color online) Schematic views of (a) Si/Si02 and (b) Ge/Ge02 models with an oxygen 
vacancy after geometrical optimization. Dashed lines indicate the boundaries of supercells. The 
key to the symbols is the same as in Fig. 



FIG. 4. (Color online) Conductance of Si/Si02 and Ge/Ge02 models as functions of incident 
electron energy measured from the valence band maximum (VBM) of the substrates. Red circles 
and black squares represent the conductance spectra of Si/Si 02 and Ge/Ge 02 models, respectively. 
VBM and conduction band minimum (CBM) of Si (Ge) substrate are indicated by vertical dashed 
(dotted) lines. 
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FIG. 5. (Color online) Conductance of Si/Si02 and Ge/Ge02 models as functions of incident 
electron energy measured from the valence band maximum (VBM) of the substrates. The key to 
the symbols is the same as in Fig. VBM and conduction band minimum (CBM) of Si (Ge) 
substrate are indicated by vertical dashed (dotted) lines. 



FIG. 6. (Color online) Contour plot of the charge-density distribution of electrons flowing over the 
Si/Si02 model with incident energy of VBM -|- 0.41 eV. The key to the symbols is the same as in 
Fig. j^a). Here the charge density is integrated in the [100] {x) direction. Each contour represents 
twice or half the density of adjacent contours; the lowest contour is 1.97 x 10“"^ electron/eV/A^. 
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FIG. 7. (Color online) Contour plots of the density of states (DOS) (left panels) and charge density 
of dangling-bond (DB) states (right panels) for (a) Si/Si02 and (b) Ge/Ge02 models without 
electrodes. Arrows denote the DOS peaks derived from the DB states. Energies are measured 
from the Fermi level Ep and the z coordinate of the atom with DB is set to zero. In the right 
panels, the charge density is integrated in the [100] (x) direction. The key to the symbols in the 
right panels is the same as in Fig. The coordinate in the z direction corresponds to that of 
the left panel. Each contour represents twice or half the density of adjacent contours; the lowest 
contour is 2.56 x 10“^ electron/eV/A (2.79 x 10“^ electron/A^) in the left (right) panels. 
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